Unsupervised learning on illicit trade data
Alice Lépissier
Bren School of Environmental Science and Management
Policy and media attention on illicit financial flows (IFF) has increased, with the recognition that Africa is a net creditor to the world.
In [3]:
listOfImageNames = ['Images/NYT_op-ed_2013.PNG',
'Images/Guardian headline_2017.PNG',
'Images/Guardian headline_2015.PNG',
'Images/Economist headline_2019.PNG']
for imageName in listOfImageNames:
display(Image(filename=imageName))
What is trade mis-invoicing?
- The deliberate mis-statement of price or quantity of internationally traded goods in invoices presented to customs
- Can occur at import or export
- Can result in an inflow or outflow of money
Motivations for trade mis-invoicing include:
- Evading tariffs
- Exploiting subsidy regimes
- Hiding transfers of wealth
Why does trade mis-invoicing matter?
- Outflows undermine the fiscal base and public spending
- Financing gap needed to meet the Sustainable Development Goals (SDGs)
- Combating trade mis-invoicing is crucial for the mobilization of domestic resources in the continent, and can catalyze sustainable development
In [4]:
display(Image(filename='Images/Governance_loop.png'))
Data source
- trade mis-invoicing data-set of the United Nations Economic Commission for Africa
- panel with $n=6,248,254$ of mis-invoiced trade between 179 reporting jurisdictions and 179 partner countries for years 2000-2016, disaggregated commodities (HS2 code)
Zoom in on Africa
- During 2000-2016, the continent lost on average \$83 billion a year in gross illicit outflows
- Net cumulative flows during that period were \$362 billion
In [5]:
display(Image(filename='Images/GER Total Average Africa IFF Percent Trade high.png'))
Goals
- Perform unsupervised learning on this data to extract meaningful insights
- Dimension reduction (PCA) - can I recover policy-relevant dimensions of variation?
Data wrangling¶
- Data-set of mis-invoiced trade by sector (unit of observation is reporter-year)
- View 1: unit of observation = reporter, features = sectors
- View 2: unit of observation = sectors, features = reporters
- Bilateral trade matrix (dyadic data)
Mis-invoiced trade for countries by sectors
In [9]:
IFF_Sector_Imp = IFF_Sector[['section', 'Imp_IFF_hi']]
IFF_Sector_Exp = IFF_Sector[['section', 'Exp_IFF_hi']]
IFF_Sector_Imp
Out[9]:
Mis-invoiced trade for dyads
In [12]:
IFF_Dest_Imp = IFF_Dest[['partner.ISO', 'Imp_IFF_hi']]
IFF_Dest_Exp = IFF_Dest[['partner.ISO', 'Exp_IFF_hi']]
IFF_Dest_Imp
Out[12]:
More data wrangling
In [14]:
crosswalk = pd.read_excel("Data/crosswalk.xlsx").rename(columns={'Country': 'country'})
crosswalk.head()
Out[14]:
Auxiliary functions¶
In [15]:
def create_features(data, values, features, obs):
"""
TODO: add doc string
"""
features_data = data.pivot_table(values=values,
columns=features,
index=[obs, 'year'],
fill_value=0)
return features_data
In [17]:
def biplot_PCA(features_data, nPC=2, firstPC=1, secondPC=2, obs='reporter.ISO', show_loadings=False):
"""Projects the data in the 2-dimensional space spanned by 2 principal components
chosen by the user, along with a bi-plot of the top 3 loadings per PC, and colors
by class label.
Args:
features_data: data-set of features
nPC: number of principal components
firstPC: integer denoting first principal component to plot in bi-plot
secondPC: integer denoting second principal component to plot in bi-plot
obs: string denoting index of class labels (in features_data)
show_loadings: Boolean indicating whether PCA loadings should be displayed
Returns:
plot (static)
pca_loadings (if show_loadings=True)
"""
# Run PCA (standardize data beforehand)
features_data_std = StandardScaler().fit_transform(features_data)
# features_data_std = power_transform(features_data, method='yeo-johnson', standardize=True)
pca = PCA(n_components=nPC, random_state=234)
princ_comp = pca.fit_transform(features_data_std)
# Extract PCA loadings
cols = ['PC' + str(c+1) for c in np.arange(nPC)]
pca_loadings = pd.DataFrame(pca.components_.T,
columns=cols,
index=list(features_data.columns))
# Extract PCA scores
pca_scores = pd.DataFrame(princ_comp,
columns=cols)
pca_scores[obs] = features_data.reset_index()[obs].values.tolist()
pca_scores['year'] = features_data.reset_index()['year'].values.tolist()
score_PC1 = princ_comp[:,firstPC-1]
score_PC2 = princ_comp[:,secondPC-1]
# Generate plot data
if obs == 'reporter.ISO':
plot_data = pd.merge(pca_scores, obs_info, on=[obs, 'year'])
color_obs = 'reporter'
tooltip_obs = ['reporter', 'year', 'Income group (World Bank)', 'Country status (UN)']
else:
plot_data = pca_scores
color_obs = 'section'
tooltip_obs = ['section', 'year']
# Return chosen PCs to plot
PC1 = 'PC'+str(firstPC)
PC2 = 'PC'+str(secondPC)
# Extract top loadings (in absolute value)
# TO DO: use dict to iterate over
toploadings_PC1 = pca_loadings.apply(lambda x: abs(x)).sort_values(by=PC1).tail(3)[[PC1, PC2]]
toploadings_PC2 = pca_loadings.apply(lambda x: abs(x)).sort_values(by=PC2).tail(3)[[PC1, PC2]]
originsPC1 = pd.DataFrame({'index':toploadings_PC1.index.tolist(),
PC1: np.zeros(3),
PC2: np.zeros(3)})
originsPC2 = pd.DataFrame({'index':toploadings_PC2.index.tolist(),
PC1: np.zeros(3),
PC2: np.zeros(3)})
toploadings_PC1 = pd.concat([toploadings_PC1.reset_index(), originsPC1], axis=0)
toploadings_PC2 = pd.concat([toploadings_PC2.reset_index(), originsPC2], axis=0)
toploadings_PC1[PC1] = toploadings_PC1[PC1]*max(score_PC1)*1.5
toploadings_PC1[PC2] = toploadings_PC1[PC2]*max(score_PC2)*1.5
toploadings_PC2[PC1] = toploadings_PC2[PC1]*max(score_PC1)*1.5
toploadings_PC2[PC2] = toploadings_PC2[PC2]*max(score_PC2)*1.5
# Project top 3 loadings over the space spanned by 2 principal components
lines = alt.Chart().mark_line().encode()
for color, i, dataset in zip(['#440154FF', '#21908CFF'], [0,1], [toploadings_PC1, toploadings_PC2]):
lines[i] = alt.Chart(dataset).mark_line(color=color).encode(
x= PC1 +':Q',
y= PC2 +':Q',
detail='index'
).properties(
width=500,
height=500
)
# Add labels to the loadings
text=alt.Chart().mark_text().encode()
for color, i, dataset in zip(['#440154FF', '#21908CFF'], [0, 1], [toploadings_PC1[0:3], toploadings_PC2[0:3]]):
text[i] = alt.Chart(dataset).mark_text(
align='left',
baseline='bottom',
color=color
).encode(
x= PC1 +':Q',
y= PC2 +':Q',
text='index'
)
# Scatter plot colored by observation class label
points = alt.Chart(plot_data).mark_circle(size=60).encode(
x=alt.X(PC1, axis=alt.Axis(title='Principal Component ' + str(firstPC))),
y=alt.X(PC2, axis=alt.Axis(title='Principal Component ' + str(secondPC))),
color=alt.Color(color_obs, scale=alt.Scale(scheme='category20b'),
legend=alt.Legend(orient='right')),
tooltip=tooltip_obs
).interactive()
# Bind it all together
chart = (points + lines[0] + lines[1] + text[0] + text[1])
chart.display()
if show_loadings:
return pca_loadings
In [18]:
def scree_plot(features_data, show_explained_var=False):
"""
TODO: add comments and docstring
"""
features_data_std = StandardScaler().fit_transform(features_data)
pca = PCA(n_components=features_data_std.shape[1], random_state=234)
princ_comp = pca.fit_transform(features_data_std)
explained_var = pd.DataFrame({'PC': np.arange(1,features_data_std.shape[1]+1),
'var': pca.explained_variance_ratio_,
'cumvar': np.cumsum(pca.explained_variance_ratio_)})
# Adapted from https://altair-viz.github.io/gallery/multiline_tooltip.html
# Create a selection that chooses the nearest point & selects based on x-value
nearest = alt.selection(type='single', nearest=True, on='mouseover',
fields=['PC'], empty='none')
# The basic line
line = alt.Chart(explained_var).mark_line(interpolate='basis', color='#FDE725FF').encode(
alt.X('PC:Q',
scale=alt.Scale(domain=(1, len(explained_var))),
axis=alt.Axis(title='Principal Component')
),
alt.Y('cumvar:Q',
scale=alt.Scale(domain=(min(explained_var['cumvar']), 1)),
axis=alt.Axis(title='Cumulative Variance Explained')
),
)
# Transparent selectors across the chart. This is what tells us
# the x-value of the cursor
selectors = alt.Chart(explained_var).mark_point().encode(
x='PC:Q',
opacity=alt.value(0),
).add_selection(
nearest
)
# Draw points on the line, and highlight based on selection
points = line.mark_point().encode(
opacity=alt.condition(nearest, alt.value(1), alt.value(0))
)
# Draw text labels near the points, and highlight based on selection
text = line.mark_text(align='left', dx=5, dy=-5).encode(
text=alt.condition(nearest, 'cumvar:Q', alt.value(' '))
)
# Draw a rule at the location of the selection
rules = alt.Chart(explained_var).mark_rule(color='gray').encode(
x='PC:Q',
).transform_filter(
nearest
)
# Put the five layers into a chart and bind the data
out = alt.layer(
line, selectors, points, rules, text
).properties(
title='Cumulative scree plot',
width=600, height=300
)
out.display()
if show_explained_var:
return explained_var[['PC', 'var']]
PCA on feature space (for individual reporting countries)¶
Sector features¶
In [19]:
sector_features = create_features(IFF_Sector_Imp, 'Imp_IFF_hi',
features='section', obs='reporter.ISO')
sector_features
Out[19]:
Biplots¶
In [20]:
biplot_PCA(sector_features, 10, 1, 2, obs='reporter.ISO', show_loadings=True)
Out[20]:
In [21]:
display(Image(filename='Images/Treemap sectors average Africa.png'))
In [22]:
biplot_PCA(sector_features, 10, 3, 4, obs='reporter.ISO')
In [23]:
biplot_PCA(sector_features, 10, 5, 6, obs='reporter.ISO')
Explained variance¶
In [24]:
scree_plot(sector_features, show_explained_var=True)
Out[24]:
Country features¶
In [31]:
country_features = create_features(IFF_Sector_Imp, 'Imp_IFF_hi',
features='reporter.ISO', obs='section')
country_features
Out[31]:
Biplots¶
In [32]:
biplot_PCA(country_features, 10, 1, 2, obs='section', show_loadings=True)
Out[32]:
In [33]:
biplot_PCA(country_features, 10, 3, 4, obs='section')
In [34]:
biplot_PCA(country_features, 10, 5, 6, obs='section')
Explained variance¶
In [35]:
scree_plot(country_features)
Variance-stabilizing and normalizing transformations¶
In [38]:
fig, axes = joypy.joyplot(country_features.iloc[:,0:22], colormap=plt.cm.magma, figsize=(8,8));
fig, axes = joypy.joyplot(country_features.iloc[:,22:47], colormap=plt.cm.magma, figsize=(8,8));
In [39]:
country_features_log = country_features.apply(lambda x: np.log(x+1) if np.issubdtype(x.dtype, np.number) else x, axis=0)
In [40]:
fig, axes = joypy.joyplot(country_features_log.iloc[:,0:22], colormap=plt.cm.magma, figsize=(8,8));
fig, axes = joypy.joyplot(country_features_log.iloc[:,22:47], colormap=plt.cm.magma, figsize=(8,8));
In [41]:
country_features_yeo = power_transform(country_features, method='yeo-johnson', standardize=True)
country_features_yeo = pd.DataFrame(country_features_yeo,
index=country_features.index,
columns=country_features.columns)
fig, axes = joypy.joyplot(country_features_yeo.iloc[:,0:22], colormap=plt.cm.magma, figsize=(8,8));
fig, axes = joypy.joyplot(country_features_yeo.iloc[:,22:47], colormap=plt.cm.magma, figsize=(8,8));
In [42]:
biplot_PCA(country_features_log, 10, 1, 2, obs='section')
In [43]:
biplot_PCA(country_features_yeo, 10, 1, 2, obs='section')
PCA on bilateral trade matrix¶
In [44]:
partner_features = create_features(IFF_Dest_Imp, 'Imp_IFF_hi',
features='partner.ISO', obs='reporter.ISO')
partner_features
Out[44]:
Biplots¶
In [45]:
biplot_PCA(partner_features, 10, 1, 2, show_loadings=True)
Out[45]:
In [46]:
display(Image(filename='Images/Top 10 destinations GER average Africa.png'))
In [47]:
biplot_PCA(partner_features, 10, 3, 4)
In [48]:
biplot_PCA(partner_features, 10, 5, 6)
Explained variance¶
In [49]:
scree_plot(partner_features)
Intra-African illicit financial flows¶
In [59]:
partner_features_AFR = create_features(IFF_Dest_Imp_AFR, 'Imp_IFF_hi',
features='partner.ISO', obs='reporter.ISO')
partner_features_AFR
Out[59]:
In [63]:
biplot_PCA(partner_features_AFR, partner_features_AFR.shape[1], 1, 2)
Remove outliers¶
In [64]:
features_data_std = StandardScaler().fit_transform(partner_features_AFR)
pca = PCA(n_components=partner_features_AFR.shape[1])
princ_comp = pca.fit_transform(features_data_std)
#outlying = (princ_comp[:,44] > 1.5) | (princ_comp[:,45] > 3e-15)
#outlying = (princ_comp[:,44] > 0.5) | (princ_comp[:,45] > 1.5e-15)
outlying = (princ_comp[:,0] > 6) | (princ_comp[:,1] > 4)
partner_features_AFR[outlying].index
Out[64]:
In [65]:
partner_features_AFR_noout = partner_features_AFR[~outlying]
In [66]:
biplot_PCA(partner_features_AFR_noout, 46, 1, 2)
Coloring by class label¶
In [69]:
biplot_PCA_classes(sector_features, 21, 1, 2, classes='UN_Intermediate_Region')
In [70]:
biplot_PCA_classes(sector_features, 21, 1, 2, classes='Income group (World Bank)')
In [71]:
biplot_PCA_classes(partner_features, 46, 1, 2)
In [72]:
biplot_PCA_classes(partner_features, 46, 1, 2, classes='Income group (World Bank)')
In [73]:
biplot_PCA_classes(partner_features_AFR, 46, 1, 2)
In [74]:
biplot_PCA_classes(partner_features_AFR, 46, 1, 2, classes='Income group (World Bank)')
In [75]:
biplot_PCA_classes(partner_features_AFR_noout, 46, 1, 2)
In [76]:
biplot_PCA_classes(partner_features_AFR_noout, 46, 1, 2, classes='Income group (World Bank)')
Heatmap¶
In [77]:
# Change column names because Vega complains when there is a period
data = IFF_Dest_AFR.reset_index().groupby(['reporter.ISO', 'partner.ISO']).mean() \
.reset_index().rename(columns={'reporter.ISO': 'r', 'partner.ISO': 'p'})\
.apply(lambda x: np.log10(x+1) if np.issubdtype(x.dtype, np.number) else x, axis=0)
alt.Chart(
data,
title='Mis-invoiced trade'
).mark_rect().encode(
alt.X('r:N', axis=alt.Axis(title='Reporting country')),
alt.Y('p:N', axis=alt.Axis(title='Partner country')),
color=alt.Color('Imp_IFF_hi:Q',
scale=alt.Scale(scheme='inferno', reverse=True, nice=True)),
tooltip=[
alt.Tooltip('Imp_IFF_hi:Q', title='Mis-invoiced imports'),
alt.Tooltip('r', title='Reporting country'),
alt.Tooltip('p', title='Partner country')
]
).properties(width=550, height=550)
Out[77]:
Concluding thoughts
- Target 16.4 of the SDGs: "By 2030, significantly reduce illicit financial and arms flows, strengthen the recovery and return of stolen assets and combat all forms of organized crime"
- Increasing international policy attention to the problem has lead to demands for more accurate estimates.
- Accurate estimates of the nature and of the scale of the problem can support policy prioritization and targeted intervention to combat illicit flows.
Next steps
- hierarchical clustering (and reordered heatmap)
- explore clustering a bit more
- NMF on bilateral trade matrix which also includes sectors, so unit of observation is dyad, and features are sectors
- Merge in GDP and (maybe?) distance measures to see if NMF recovers gravity interpretation of international trade